%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AGUtmpl.tex: this template file is for articles formatted with LaTeX2e,
% Modified July 2011
%
% This template includes commands and instructions
% given in the order necessary to produce a final output that will
% satisfy AGU requirements.
%
% PLEASE DO NOT USE YOUR OWN MACROS
% DO NOT USE \newcommand, \defcommand, or \renewcommand.
%
% FOR FIGURES, DO NOT USE \psfrag or \subfigure.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% All questions should be e-mailed to latex@agu.org.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Step 1: Set the \documentclass
%
% There are two options for article format: two column (default)
% and draft.
%
% PLEASE USE THE DRAFT OPTION TO SUBMIT YOUR PAPERS.
% The draft option produces double spaced output.
%
% Choose the journal abbreviation for the journal you are
% submitting to:

% jgrga JOURNAL OF GEOPHYSICAL RESEARCH
% gbc   GLOBAL BIOCHEMICAL CYCLES
% grl   GEOPHYSICAL RESEARCH LETTERS
% pal   PALEOCEANOGRAPHY
% ras   RADIO SCIENCE
% rog   REVIEWS OF GEOPHYSICS
% tec   TECTONICS
% wrr   WATER RESOURCES RESEARCH
% gc    GEOCHEMISTRY, GEOPHYSICS, GEOSYSTEMS
% sw    SPACE WEATHER

% (If you are submitting to a journal other than jgrga,
% substitute the initials of the journal for "jgrga" below.)

\documentclass[draft,jgrga]{AGUTeX}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OPTIONAL:
% To produce a two-columned version:
% \documentclass[jgrga]{AGUTeX}

% Two-columned format can be used to estimate the number of pages
% for the final published PDF.

% PLEASE USE THE DRAFT OPTION TO SUBMIT YOUR PAPERS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OPTIONAL:
% To create numbered lines:

% If you don't already have lineno.sty, you can download it from
% http://www.ctan.org/tex-archive/macros/latex/contrib/ednotes/
% (or search the internet for lineno.sty ctan), available at TeX Archive Network (CTAN).
% Take care that you always use the latest version.

% To activate the commands, uncomment \usepackage{lineno}
% and \linenumbers*[1]command, below:

% \usepackage{lineno}
% \linenumbers*[1]

%  To add line numbers to lines with equations:

%  \begin{linenomath*}
%  \begin{equation}
%  \end{equation}
%  \end{linenomath*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures and Tables
%
% When submitting articles through the GEMS system:
% COMMENT OUT ANY COMMANDS THAT INCLUDE GRAPHICS.
% (See FIGURES section near the end of the file.)
%
% DO NOT USE \psfrag or \subfigure commands.
%
%  Figures and tables should be placed AT THE END OF THE ARTICLE,
%  after the references.
%
%  Uncomment the following command to include .eps files
%  (comment out this line for draft format):
%  \usepackage[dvips]{graphicx}
%
%  Uncomment the following command to allow illustrations to print
%   when using Draft:
%  \setkeys{Gin}{draft=false}
%
% Substitute one of the following for [dvips] above
% if you are using a different driver program and want to
% proof your illustrations on your machine:
%
% [xdvi], [dvipdf], [dvipsone], [dviwindo], [emtex], [dviwin],
% [pctexps],  [pctexwin],  [pctexhp],  [pctex32], [truetex], [tcidvi],
% [oztex], [textures]
%
% See how to enter figures and tables at the end of the article, after
% references.
%
% Use units package to manage the descriptions of units in the text.
% This is non-standard for the AGU template, and added by JVJ.
\usepackage[thickqspace,mediumspace,squaren]{SIunits}
%% ------------------------------------------------------------------------ %%
%
%  ENTER PREAMBLE
%
%% ------------------------------------------------------------------------ %%

% Author names in capital letters:
\authorrunninghead{SACK AND JOHNSON}

% Shorter version of title entered in capital letters:
\titlerunninghead{Propagation of errors in bed elevation}

% Author mailing address: please repeat this command for
% each author and alphabetize authors:

\authoraddr{K. Sack,
Department of Civil Engineering, CERECAM University of Cape Town, ZA.
(kev.sack@gmail.com)}

\authoraddr{J. V. Johnson,
Department of Computer Science, University of
Montana, 32 Campus Ave., Missoula, MT 59801, USA.
(jesse.v.johnson@gmail.com)}

\begin{document}

%% ------------------------------------------------------------------------ %%
%
%  TITLE
%
%% ------------------------------------------------------------------------ %%


\title{The propagation of errors in bed elevation to adjoint based assimilation
of surface velocity}
%
% e.g., \title{Terrestrial ring current:
% Origin, formation, and decay $\alpha\beta\Gamma\Delta$}
%

%% ------------------------------------------------------------------------ %%
%
%  AUTHORS AND AFFILIATIONS
%
%% ------------------------------------------------------------------------ %%


%Use \author{\altaffilmark{}} and \altaffiltext{}

% \altaffilmark will produce footnote;
% matching \altaffiltext will appear at bottom of page.

\authors{K. Sack,\altaffilmark{1}
J. V. Johnson,\altaffilmark{2}}

\altaffiltext{1}{Department of Mechanical Engineering, CERECAM University of Cape Town}

\altaffiltext{2}{Department of Computer Science, University of Montana}


%% ------------------------------------------------------------------------ %%
%
%  ABSTRACT
%
%% ------------------------------------------------------------------------ %%

% >> Do NOT include any \begin...\end commands within
% >> the body of the abstract.

\begin{abstract}
This work is presented in two parts. In the first, we use a numerical inversion
technique to determine a ``mass conserving bed'' (mcb). This technique adjusts the bed
elevation to assure that the mass flux determined from surface velocity
measurements does not violate conservation. In this way, the mcb becomes an
interpolation algorithm, allowing us to estimate the bed between radar flight
lines. We then perform cross validation on this technique by using a subset of
available flight lines; the unused flight lines provide data to compare to, and
quantify the errors produced by the mcb interpolation method. We find that the
errors are similar to those produced with more conventional interpolation
schemes, such as kriging, but in this case the bed interpolation is much more
consistent with the physics that govern ice sheet models. 

In the second part, a numerical model including the full, non-Newtonian
momentum balance of glacial ice used to relate the errors in basal geometry to
the kinematic surface boundary condition via inversion of surface velocities to
find the basal traction. The model domain is a two-dimensional, periodic,
flow-line that includes an obstacle at the bed. A control run is carried out to
establish the surface velocity produced by the model. The control surface
velocity is then used as a target for data inversions performed on perturbed
versions of the original control beds. In this way, errors in bed measurement
are propagated forward to investigate errors in the evolution of the free
surface. The growth of errors with perturbation size is demonstrated to be
exponential for cases in which the perturbation is positive, or producing a
hill. In cases where the perturbation is negative, or forming a depression in
the domain, errors grow exponentially, but then decrease. We argue that this is
due to the depression being deep enough ($\sim$ 80 m) that ice shears over the
top of the depression, effectively returning the bed to something close to the
unperturbed state. Our primary conclusion relates the size of errors in the
surface evolution to errors in the bed, and suggests an optimal spacing of
flight lines.  \end{abstract}

%% ------------------------------------------------------------------------ %%
%
%  BEGIN ARTICLE
%
%% ------------------------------------------------------------------------ %%

% The body of the article must start with a \begin{article} command
%
% \end{article} must follow the references section, before the figures
%  and tables.

\begin{article}

\section{Introduction}

One way the consequences of anthropogenic sea-level rise \citep{Solomon:2007} is
being assessed is by collection of three dimensional geometry information for
the Antarctic and Greenland ice sheets. The largest such effort is NASA's
Operation IceBridge \cite{Studinger:2010bl}, which is leading to the
development of gridded, three-dimensional geometric information at the 500-1000
meter level \citep{SeaRISE}. While enticing, geometry alone does not provide the
insight into ice dynamics required for assessment. Further information is
provided by surface velocity data derived from Interferometric Synthetic
Aperture Radar (InSAR) \citep[e.g.][]{Joughin:2010th}.  Existing on grids similar or
superior in resolution to the geometric data, surface velocity data reveal
complex networks working to drain accumulation deposited on ice sheets \citep{Rignot:2011}.

In order to extend the spatial and temporal coverage of datasets, many look to
numerical models. If models are to reliably estimate the relation between 
warmer climate and sea-level, they will require detailed geometric and boundary
condition data. How detailed data a major concern, because the cost of any data
collection effort is significant, and there exists a constant tension between
depth vs.  breadth of coverage. There is likely some optimum for predictive
modeling; a point where additional data do not significantly alter the modeled
flow. In the case of the geometric information being acquired by Operation
IceBridge, locating this optimum first requires assessment of errors in geometry
that can be tolerated by numerical models. After that, the errors associated
with interpolation between flight lines can be assessed and the appropriate
spacing of flight lines be determined.

Hence, in this article we explore the significance of errors in the bed
elevation with respect to an popular means of initializing ice sheet models;
assimilation of surface velocity data \citep{Larour:2012fk,Brinkerhoff:2011vr}.
This assimilation process alters the basal traction in order to recover
velocities observed on the surface \citep{MacAyeal:1992a}. Because the basal
traction field can take arbitrary values at each point in the domain, the
optimization procedure reproduces surface velocity data to within narrow
tolerances \citep{Morlighem:2010cr}. However, questions remain about how closely
the inverted traction field represents a physical quantity.  For example, the
optimized basal traction field may be compensating for errors in the observed
geometry, introducing errors in the traction field and misrepresenting the
physics that dominate glacier sliding. Such errors are problematic because a
system initialized via assimilation will evolves in time, moving away from an
inconsistent initial state, and relaxing to one that is consistent with the
model's physics. This process introduces artificial transients that will
strongly bias model results. 

Here, the severity of problems related to assimilating data onto imperfect
geometry is assessed with a set of numerical modeling experiments. Our strategy
is to establish a control run for comparison to a set of runs having perturbed
geometries. All runs on perturbed geometries involve assimilation of the surface
velocity from the control run, emulating what would happen with real data
efforts which use relatively high accuracy surface velocity data but imperfect
geometries \cite{Seroussi:2011vy}.

The assimilation methods used here were pioneered by \citep{MacAyeal:1992a}, and
seen considerable use as more surface velocity data has become available,
\citet[e.g.][]{Rignot:2011}.  \citet{Seroussi:2011vy} were the first to document
how errors in the bed geometry can play a significant role in the assimilation
of data, and lead to problems with the conservation of mass. Other
investigations of the errors associated with the bed are found in XXX and XXX.
This work is the first to systematically probe the effects of errors in the bed
on surface evolution.

We use numerical modeling, as the relations are too complex to be tractable in
closed form analysis. More words on this...


\section{Description of Experiments}

For the sake of comparison, we assume perfect knowledge of the geometry,
physics, and basal traction for some control configuration. When the model is
run on the control configuration, its output will become the 'data' that
subsequent runs are forced to match via optimization. Subsequent model runs are
performed on geometries that are not the same as the control; their beds are
systematically perturbed at some point. Yet, using data assimilation, control
run surface velocities are reproduced to within some tolerance. 

At this point, the model has done exactly what was desired; errors in bedrock
geometry are now part of the basal traction field. At this point, a metric for
assessing the significance of the errors is needed, and we opt for the changes
in the kinematic boundary condition,

\begin{equation}
\frac{\partial s} {\partial t} = -\mathbf{u} \cdot \nabla s + a.
\label{eq:KBC}
\end{equation}

This boundary condition dictates the time
evolution of the surface of the ice sheet, and errors in it capture how errors
in the bed propagate forward in time. Again, we pretend we know what this should
be, based on output from the control run. Hence, results are reported as a sum
squared error difference with the control run, as well as plots on the same axes
to show the spatial structure of errors.

\subsection{Model Geometry}
To create a realistic bed topography we apply a numerical algorithm to synthesize
realistic bed roughnesses. The algorithm is only applied in the center of the bed
topography to preserve and smooth the geometry leading to and emanating from our
periodic boundary conditions (as depicted in fig \ref{fig:setup2}).\\

Given the domain of the terrain [0,L], the end points of the function $f(0)$ and
$f(L)$, a perturbation maximum $h$ and a regularization parameter $r$ we calculate
$f(mid)$ using a Fractal Brownian Motion (FBM) Method \citep{Fournier}. Once $f(mid)$ is
calculated we subdivide between known points and calculate new \textit{midpoint}
values between points. \\

%\begin{figure}[h]
%\centering
%		\includegraphics[scale=0.35]{pictures/3stepa.png}
%	\caption{Basic results of FBM method. The number of known points after $n$ iterations is $2^n +1$.}
%	\label{fig:mid3}
%\end{figure}

A key principle to the FBM is that with each further subdivision the random
displacements should decrease by a factor related to the prescribed regularization 
parameter $r$, which controls roughness. The most
commonly used factor is given by $2^{-r}$, and $r$ is limited to the range
$[0,1]$. As $r \rightarrow 1$, the factor $2^{-r} \rightarrow 1/2$ which
corresponds to halving the potential displacement in each iteration.\\

The following result below in (fig.\ref{fig:Roughness}) tracks closely with the literature [XXX].\\

\begin{figure}[h]
	\centering
		\includegraphics[width=0.45\textwidth]{pictures/Roughness.png}
	\caption{Basic results of Midpoint displacement method for 8 subdivisions with varying roughness factor}
	\label{fig:Roughness}
\end{figure}

Once the fractal terrain was generated, it was imported as a set of
displaced values as bed geometry into the model. To ensure the domain is atleast $C^1$ continuous, the fractal bed data set is interpolated using a series of \emph{b\`{e}zier} cubic spline functions. One can consider the entire bed generation process  summarised by figure \ref{fig:Bedcreation}.

\begin{figure}[h]
	\centering
		\includegraphics[scale=0.4]{AGU_RES/bedcreation.png}
	\caption{Process invovled in creating the set of synthetic bed data. The upper figures correspond to to the fractal generation routine and the bottom figure corresponds to the $C^1$ continuous interpolated bed}
\label{fig:Bedcreation}
\end{figure}


All finite element modelling is done in the \emph{COMSOL Multiphysics} environment, a commercial software package. The terrain generation modelling was created using \emph{MATLAB}.\\

%The following figure is the idealized reference geometry we will use to base our experiments on:
%
%\begin{figure}[h]
%	\centering
%		\includegraphics[scale = 0.4]{graphics/mesh.png}
%	\caption{Discretized 2D geometry with 756 nodes}
%	\label{fig:2dmesh}
%\end{figure}

\subsection{Perturbations of geometry}

In order to investigate the model's ability to assimilate data using incorrect parameters, we introduce perturbations into the bed geometry. We limit the scope of our investigation to single perturbations in each experiment, increasing the perturbation size from $0 - 150m$ in each successive experiment. As the bed is created using a cubic spline functions, a single perturbation at point will affect the localized domain surrounding the point for a radius between $2-3km$, resulting in an evenly distributed and smooth perturbation to an area of the geometry as opposed to a sharp perturbation of a single point. 

\begin{figure}[h]
	\centering
		\includegraphics[scale = 0.3]{AGU_RES/bedperturbed.png}
	\caption{Illustration of a downwards perturbation and the influence to the surrounding beg geometry.}
	\label{fig:perturbed}
\end{figure}


We explore both negative (ie. valley forming) and positive (ie. hill forming) perturbations.

\section{Description of the Numerical Model}

\begin{table}

\begin{center}

\begin{tabular}{l|l|l}

\textbf{Parameter}& \textbf{Symbol} &\textbf{Value}\\
\hline
\hline
Gravitational acceleration&\textbf{g}&9.81 \metre\per\second\squared\\

Density of ice & $\rho$ & 911 \kilogrampercubicmetre\\

Seconds per year & -& 31556926 \second\usk a$^{-1}$\\

Glen's flow law exponent & $n$ & 3\\

Viscosity regularization & $\dot \epsilon$ & 10$^\mathrm{-30}$ \pascalsecond \\
\hline
\end{tabular}

\end{center}
\caption{Parameters and physical constants used in the model.}

\end{table}

\subsection{Field equations}
Our model is built upon the continuum mechanical formulation of the laws of
conservation of mass and momentum for an incompressible fluid. These
are, respectively;

\begin{equation}
\label{incompres}
\nabla \cdot \mathbf{u} = 0,
\end{equation}
\begin{equation}
\label{momentum}
\rho \frac{d\mathbf{u}}{dt} = \nabla \cdot \sigma + \rho \mathbf{g},
\end{equation}

$\mathbf{u}$ represents the velocity vector, and $\sigma$ the stress tensor.
Physical constant $\mathbf{g}$ is defined in Table ~\ref{params}.  Analysis is
restricted to the $xz$ plane, or the vertical profile, making $\nabla =
\frac{\partial}{\partial x}\hat{i} + \frac{\partial}{\partial z}\hat{k}$, where
$\hat{i}$ and $\hat{k}$ are unit vectors in the $x$ and $z$ directions,
respectively.  

\subsubsection{Conservation of momentum and mass}
The constitutive relation for ice takes the form
\begin{equation}
\tau_{ij} = 2 \eta \dot{\epsilon}_{ij},
\end{equation}
$\tau_{ij}$ is the $ij$ element of the deviatoric stress tensor, which is
defined by $\tau_{ij}=\sigma_{ij}-p\delta_{ij}$, and $\delta_{ij}$ is the
Kronecker delta function. Isotropic pressure is defined as
$p=-\frac{1}{3}\sum_i \sigma_{ii}$.  ~$\dot{\epsilon}_{ij}$ represents the
corresponding element of the strain rate tensor, and $\eta$ the viscosity. 
The strain rate tensor is given by, and related to velocity gradients as follows
\begin{equation}
\dot{\epsilon}_{ij} = \frac{1}{2} \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right).
\end{equation}

A non-Newtonian rheology is used for ice
\begin{equation}
\label{viscosity}
\eta = \frac{1}{2} A(\theta^*) ^{-1/n} (\dot{\epsilon}_\Pi +
\dot{\epsilon}_0)^{(1-n)/n},
\end{equation}
with $\dot{\epsilon}^2_{\Pi} =
\frac{1}{2}\Sigma_{ij}\dot{\epsilon}_{ij}\dot{\epsilon}_{ij} $, or the second
invariant of the strain rate tensor, and $\dot{\epsilon}_0$ is a regularization
parameter introduced to avoid a singularity at zero strain rate.  Glen's flow
law \citep{Paterson:1994} gives $n$=3.
$A(\theta^*)$ is the flow law rate factor, given by \citet{Paterson:1994}.
where $\theta^*$ is the homologous temperature, defined by $\theta^* = \theta +\beta p$. In this work, we disregard the temperature dependence of flow and assume that $A= $ 1 $\times$ 10$^{-16}$ UNITS, which corresponds to ice at about -2 degrees C (CHECK).

Under the assumption of steady state, the velocity of the ice is then determined
from Stoke's flow confined to the $xz$ plane 

\begin{equation}
\nabla \cdot \sigma = \rho \mathbf{g},
\end{equation}

and the conservation of mass $\nabla \cdot \mathbf{u} = 0$.  

\subsection{Boundary Conditions}

\subsubsection{Conservation of Momentum and Mass Boundary Conditions}
The surface of the glacier upholds the neutral or stress free boundary
condition
\begin{equation}
\sigma \hat{\mathbf{n}} = 0,
\end{equation}
where $\hat{\mathbf{n}}$ is the outward normal unit vector.

The bed of the glacier is subjected to a Weertman style sliding law, where basal
velocity and shear stress are related as
\begin{equation}
\mathbf{\tau_b} = \beta^2 \mathbf{u}
\end{equation}

where $\beta^2$ is a positive scalar, spatially variable parameter representing
the magnitude of frictional forces at the bed, $\mathbf{\tau_b}$ is given by
\begin{equation}
\mathbf{\tau_b} = \sigma \hat{\mathbf{n}} 
\end{equation}
evaluated at the base of the glacier. Here, $\mathbf{u} \cdot \hat{\mathbf{n}}
= 0$.

The vertical boundary, imposed on the left and right sides of this boundary, are
periodic

\begin{equation}
\hat{\mathbf{n}} \cdot \mathbf{u} = 0
\end{equation}
\begin{equation}
 \sigma \cdot \hat{\mathbf{t}}= 0,
\end{equation}
where $\hat{\mathbf{t}}$ is the unit vector tangent to the divide.


\subsection{Numerical Considerations}
The model uses the finite element method to solve the field equations subject to
the boundary conditions.  Lagrange quadratic elements are used
\citep{Hughes:2000}, allowing second derivatives of the velocity to be computed.
The non-linearity resulting from the viscosity (Equation~\ref{viscosity}) is
resolved by using the modified Newton's method iterative solver
\citep{Deuflhard:1974}. The resulting linear systems were solved with UMFPACK
\citep{Davis:2004}. Model specific parameters are summarized in
Table~\ref{solver}. All numerical work was carried out in the \textit{Comsol
Multiphysics} modeling environment, a commercial package for finite element
analysis of general partial differential equations. 
%-----------------------------
% Model Numerics
%-----------------------------
\begin{table}[h]

\begin{center}

\begin{tabular}{l|l}

\textbf{Quantity}&\textbf{Value}\\
\hline
\hline
Mesh Elements & 7313 \\

Degrees of Freedom & 33998\\

 Lagrange Element Type & Quadratic\\

Initial Damping Factor & 1$\times 10^{-4}$\\

Minimum Damping Factor & 1$\times 10^{-8}$ \\

%Convergence Criterion of non-Linear Solver & $<$1$\times 10^{-6}$ \\

%Convergence Criterion of Optimization Solver & $<$1$\times 10^{-4}$ \\
\hline

\end{tabular}

\end{center}

\caption{\label{solver} Quantities of importance for model numerics of a  typical study experiment.}

\end{table}

Where the Convergence Criterion of the non-Linear Solver  was set to $<$1$\times 10^{-6}$ and the 
Convergence Criterion of the Optimization Solver at & $<$1$\times 10^{-3}$.

\subsection{Modeling assumptions}
Several assumptions were made in the development of this model, and results
must be understood with these in mind.  These assumptions are as follows:

\begin{enumerate}
\item The synthetic data sets used to generate the model geometry are
sufficiently realistic to capture the response of ice sheet models to errors in
the bed.

\item Investigation of errors at a single point in the domain are suggestive of
the types of errors that would occur if the errors were isotropically
distributed across the entire domain.

\item As in any profile type model, assume that stresses acting transverse to
the dominant flow direction are small and would not alter experimental results. 

\item Solutions which result from  the data assimilation process are sufficient
for probing the sensitivities of the system with respect to changes in the
geometry. A more complete treatment would entail the evolution of the free
surface to determine the ultimate outcome of the perturbation, but that is
beyond the scope of this work.

\end{enumerate}


\subsection{Data assimilation and model initialization}

When modelling ice dynamics there are two issues that must be addressed before
numerical experiments can be conducted. First, fields which have not been
directly measured but are significant in computing flow must be estimated.
For instance, the internal distribution of temperatures are critical
to ice dynamics, but are at best known at a few boreholes. We will refer to this
process as model initialization. Secondly, the initialized model should be in
agreement with measurements that are available. We refer to this as data
assimilation.

Our strategy in this paper will be to use steady state solutions to
conservation equations to initialize the model subject to the constrains
introduced by the data assimilation process. This is not a new idea,
\citet{MacAyeal:1993a} introduced control methods in the context of ice sheet
modelling. Here, we extend the concepts to solutions which incorporate the full
flow-line stress balance.

For data assimilation, we use the adjoint of the linear operator to compute
derivatives of an objective function, and use those slopes to minimize the
function. We have defined an objective function in terms of difference
between the observed, $u^{\mathrm{obs}}(\mathbf{x})$, and modelled,
$u^{\mathrm{mod}}(\mathbf{x})$ surface velocities,
\begin{equation}
g(u,p) = \sum_{i=1}^N (u^{\textrm{obs}}(\mathbf{x_i}) - u^{\textrm{mod}}(\mathbf{x_i}))^2
\end{equation}
which will be differentiated with respect to a parameter $p$ that we vary in
order to minimize the objective function. In this case the parameter will be
$p=\beta(x)^2$ or the basal traction. Our introduction follows that of
\citet[pages 678-684]{Strang:2007uj}.\\

It was found that initial solutions to the posed optimization problem favoured high 
frequency solutions, which in turn, model rapidly changing ice flow behavior on the 
bed. To support the more realistic, gradual changes in ice flow dynamics we introduce 
a Tikhonov regularization parameter that would penalise solutions based on gradients 
and favour smoother solutions. Our usage of the Tikhonov regulization 
rephrases the optimization problem as:

\begin{equation}
g(u,p) = \sum_{i=1}^N (u^{\textrm{obs}}(\mathbf{x_i}) - u^{\textrm{mod}}(\mathbf{x_i}))^2 + \lambda \Vert ~L p \Vert ^2
\end{equation}


Where $\lambda $ is an appropiately chosen scaling parameter and $~L$ is the Tikhonov 
regularization operator. For our specific problem $~L$ is set to be a first order gradient operator. \\

`Chain rule' differentiation yields
\begin{equation}
\frac{dg}{dp} = \frac{\partial g}{\partial u}\frac{\partial u}{\partial p} + \frac{\partial g}{\partial p},
\end{equation}
where $u$ is a solution vector containing both velocities and temperatures. The
key to efficient calculation of the derivatives of the objective function is writing 
\begin{equation}
\frac{\partial g}{\partial u} = c^T
\end{equation}
or, recognizing that the objective function is linear in $u$. It is now possible to
write the gradient as 
\begin{equation}
\frac{dg}{dp} = c^T\frac{\partial u}{\partial p} + \frac{\partial g}{\partial p} = c^TA^{-1}\frac{\partial b}{\partial p} + \frac{\partial g}{\partial p},
\end{equation}
where that  $c^TA^{-1}$ is the result of solving the ``adjoint'' linear system
$A^T \lambda = c$ for $\lambda^T = c^TA^{-1}$. Note that the original problem is
assumed to be represented by the system of linear equations $Au=b$. Hence, the
gradient for each step of an optimization algorithm (we use quasi-Newton)
requires a single extra linear solve, rather than a linear solve for each of the
many parameters, $p$. This savings makes it possible to do large inverse
problems, such as computing a basal traction for each point in the model domain
(see Fig. 2). Figure 3 then corresponds to the initialized velocity and
temperature field, or the steady state solutions to the field equations that
assimilate the data. This will provide the starting point for all numerical
experiments. In some cases, such as determination of the sensitivity to $q_g$,
the entire assimilation/initialization process is repeated with different
values.

\subsection{Summary of the experimental design}
We can consider our model, $\mathcal{M} \left( B, s, \beta^2 \right) \rightarrow u_s $, as solving for the steady state problem given a prescribed Bed, $B$, surface, $s$ and traction field along the bed, $\beta^2$ and providing the user with, amongst other measures, the surface velocity, $u_s$. We investigate the errors produced when the inverse of the model is run on a perturbed bed, $B^*$, using the same surface velocity that was solved for on the forward run of the model. In this set up, the inverse model, $\mathcal{M}^{-1} \left( u_s, s, B^* \right) \rightarrow \beta^{* 2} $, produces the data assimilated traction field, $\beta^{* 2}$. \\

In this way, the model optimized the basal traction $\beta(x)^2$ to ensure the perturbed system retains the original surface velocity, which would be a data set of observed velocities, $u^{obs}$. To investigate how accurately the model captures the physical attributes of the original system, or how well changes in the traction field, $\beta(x)^2$, can compensate for errors in the bed geometry, we look at the differences in the Kinematic boundary condition, (\ref{eq:KBC})  from our data assimilated solution and the ``observed'' reference state, integrated over the surface of our ice sheet. These are measured in the following manner:

\begin{equation}
\int{\mathrm{abs}\left( {\frac{\partial s}{\partial t}}^{obs} - {\frac{\partial s}{\partial t}}^{mod} \right)} \mathrm{d}\Gamma 
\label{eq:KBCerr}
\end{equation}

The total change in (\ref{eq:KBCerr}) reflects the absolute rate of change in volume. We normalize this error metric in a more meaningful way: We take the resulting error in (\ref{eq:KBCerr}) and calculate the incorrectly introduced mass this would add to the system over a $10$ year period. We then normalize this amount as a percentage of the total mass in the ice sheet. Subsequently, this large data set of error with corresponding perturbation size and direction will form the basis of further analysis. 

\section{Results}


$\mathcal{M}^{-1} $ was run for both upward, hill forming, and downward, valley forming, perturbations in incremental size.   At each experiment, the surface profiles throughout the ice sheet (an example is given in figure \ref{fig:VelSample}), the error profile in KBC along the surface  and finally our error metric given in (\ref{eq:KBCerr}) was captured.

\begin{figure}[H]
	\includegraphics[scale=.4]{AGU_RES/0v50.png} 
	\caption{Comparison of Velocity profile produced by $\mathcal{M}^{-1}$ performed on a 50m perturbed downward experiment}
\label{fig:VelSample}
\end{figure}

The corresponding error in the KBC along the surface is presented in figure \ref{fig:KBC50}. We draw emphasis to the oscillatory nature in figure \ref{fig:KBC50} as the motivation for (\ref{eq:KBCerr}). The inclusions of the absolute operator is essential in capturing the true error in the ice sheet.\\

\begin{figure}[h]
	\includegraphics[scale=.4]{AGU_RES/KBC50.png} 
	\caption{KBC error profile along the ice sheet surface produced by $\mathcal{M}^{-1}$ performed on a 50m perturbed downward experiment}
\label{fig:KBC50}
\end{figure}

As one would expect the resulting errors in the Kinematic boundary condition produced by  $\mathcal{M}^{-1} \left( u_s, s, B^* \right) \rightarrow \beta^{* 2} $ increase with perturbation size. The error trend as a function of perturbation size can be approximated quite accurately with linear polynomials as illustrated in figure \ref{fig:Errors} and summarized in table \ref{curves}.

\begin{figure}[H]
	\includegraphics[scale=.5]{AGU_RES/Linfit2.png} 
	\caption{Percentage of absolute error as a percentage error of volume introduced into the ice sheet over a 10 year period.}
\label{fig:Errors}
\end{figure}


\begin{table} [H]

\begin{center}

\begin{tabular}{l|l|l}

\textbf{Pert. Data set}&\textbf{Downward }&\textbf{Upward}&\\

 %\textbf{Curve fitted} &\textbf{SSE} &\textbf{R-squared} &\textbf{RMSE} &\textbf{Adjusted R-sq}\\
\hline
\hline
Curve fitted & Linear & Linear\\
SSE & 0.00124 & 0.001575\\
R-squared & 0.9943 & 0.9977\\
RMSE & 0.009763 & 0.011006  \\
Adjusted R-sq & 0.9939 & 0.9975\\

% Downward perturbations & Linear & 0.00124 & 0.9943 & 0.009763 & 0.9939\\

% Upward perturbations & Linear & 0.001575 & 0.9977 & 0.011006 & 0.9975\\
\hline
\end{tabular}
\end{center}
\caption{\label{curves} Curve fitting parameters}
\end{table}


\section{Discussion}
While the linear curve fitting results do capture the overall error trend, we notice that they don't track the error behaviour for small perturbations particularly well. Additionally, in this region of small perturbations (i.e $<25m$) we can  draw attention to the way both positive and negative results are almost identical, tracking closely together.\\

We propose that for small errors introduced in the bed geometry, the physical flow of the system remains relatively unchanged. A perturbation in either direction creates the \emph{same} global error.  However once a certain threshold for errors introduced is reached, the nature of the flow is disrupted to the extent that it creates two distinctly different error results. This ``Snap through'' event is different for positive and negative perturbations. \\

The error curve for downward perturbations was significantly less than we had originally expected. This under-performance in error results from a physical response created by  $\mathcal{M}^{-1} \left( u_s, s, B^* \right)$. The depression caused by the downward perturbation, when large enough, creates an region of stagnant flow. This region of stagnant flow affectively reduces the impact of the perturbation allowing for a $150m$ downward perturbation to produce the same amount of error in ice mass as a $90m$ perturbation upwards.\\

  
  \begin{figure}[H]
	\includegraphics[scale=0.4]{AGU_RES/snapthrough.png}\\
  \caption{Illustration of a stagnant flow created in a 150m downward perturbation. The upper figure corresponds to the reference experiment. }
  \label{fig:stagnant}
  \end{figure}

To utilize the results in a meaningful way, we combine the work presented in [J V Johnson et al 2013] that links the errors introduced in bed data from raw track line data. Combining these, we can ultimately determine ideal track line spacing to be used for raw model input based on the numerical technique chosen for estimating bed elevation data and an prescribed confidence interval for the resulting model. To impose a model which produces accurate volume estimations for ice sheet, we can limit volume change and obtain a corresponding limit to acceptable perturbations. (i.e. a limit to $10\%$ in volume change corresponds to perturbations of roughly $25m$). According to  [J V Johnson et al 2013], the maximum space between data track lines one could use would be less $1.2km$ for Kriging methods and $2km$ for Mass Balance Conservation techniques.

\section{Conclusions}

The work presented shows a usable and clear relationship between perturbations and resulting error in an Ice sheet models.
etc etc

%%% End of body of article:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optional Appendix goes here
%
% \appendix resets counters and redefines section heads
% but doesn't print anything.
% After typing  \appendix
%
% \section{Here Is Appendix Title}
% will show
% Appendix A: Here Is Appendix Title
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Optional Glossary or Notation section, goes here
%
%%%%%%%%%%%%%%
% Glossary is only allowed in Reviews of Geophysics
% \section*{Glossary}
% \paragraph{Term}
% Term Definition here
%
%%%%%%%%%%%%%%
% Notation -- End each entry with a period.
% \begin{notation}
% Term & definition.\\
% Second term & second definition.\\
% \end{notation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ACKNOWLEDGMENTS

\begin{acknowledgments}
(Text here)
\end{acknowledgments}

%% ------------------------------------------------------------------------ %%
%%  REFERENCE LIST AND TEXT CITATIONS
%
% Either type in your references using
% \begin{thebibliography}{}
% \bibitem{}
% Text
% \end{thebibliography}
%
% Or,
%
% If you use BiBTeX for your references, please produce your .bbl
% file and copy the contents into your paper here.
%
% Follow these steps:
% 1. Run LaTeX on your LaTeX file.
%
% 2. Run BiBTeX on your LaTeX file.
%
% 3. Open the new .bbl file containing the reference list and
%   copy all the contents into your LaTeX file here.
%
% 4. Comment out the old \bibliographystyle and \bibliography commands.
%
% 5. Run LaTeX on your new file before submitting.
%
% AGU DOES NOT WANT a .bib or a .bbl file. Please copy in the contents of your .bbl file here.
\bibliographystyle{plain}
\bibliography{References}


%\begin{thebibliography}{}

%\bibitem[{\textit{Kilby}(2008)}]{jskilby}
%Kilby, J. S. (2008), Invention of the integrated circuit, {\it IEEE
%Trans. Electron Devices,} \textit{23}, 648--650.

%\bibitem[{\textit{Kilby et al.}(2008)}]{jskilbye}
%Kilby, J. S., S. Smith, and R. Jones (2008), Invention of the
%integrated circuit, {\it IEEE Trans. Electron Devices,} \textit{23},
%648--650.

%\end{thebibliography}

%Reference citation examples:

%...as shown by \textit{Kilby} [2008].
%...as shown by {\textit  {Lewin}} [1976], {\textit  {Carson}} [1986], {\textit  {Bartholdy and Billi}} [2002], and {\textit  {Rinaldi}} [2003].
%...has been shown [\textit{Kilby et al.}, 2008].
%...has been shown [{\textit  {Lewin}}, 1976; {\textit  {Carson}}, 1986; {\textit  {Bartholdy and Billi}}, 2002; {\textit  {Rinaldi}}, 2003].


%...as shown by \citet{jskilby}.
%...as shown by \citet{lewin76}, \citet{carson86}, \citet{bartoldy02}, and \citet{rinaldi03}.
%...has been shown \citep{jskilbye}.
%...has been shown \citep{lewin76,carson86,bartoldy02,rinaldi03}.
%
% Please use ONLY \citet and \citep for reference citations.
% DO NOT use other cite commands (e.g., \cite, \citeyear, \nocite, \citealp, etc.).

%% ------------------------------------------------------------------------ %%
%
%  END ARTICLE
%
%% ------------------------------------------------------------------------ %%

\end{article}

%% Enter Figures and Tables here:

% When submitting articles through the GEMS system:
% COMMENT OUT ANY COMMANDS THAT INCLUDE GRAPHICS.
%
% FOR FIGURES, DO NOT USE \psfrag or \subfigure commands.
%
% Figure captions go below the figure.
% Table titles go above tables; all other caption information
%  should be placed in footnotes below the table.

% DRAFT figure/table, including eps graphics
%
% \begin{figure}
% \noindent\includegraphics[width=20pc]{samplefigure.eps}
% \caption{Caption text here}
% \end{figure}
% \end{document}
%
% \begin{table}
% \caption{}
% \end{table}
%
% ---------------
% TWO-COLUMN figure/table
%
% \begin{figure*}
% \noindent\includegraphics[width=39pc]{samplefigure.eps}
% \caption{Caption text here}
% \end{figure*}
%
% \begin{table*}
% \caption{Caption text here}
% \end{table*}
%
% ---------------
% EXAMPLE TABLE
%
%\begin{table}
%\caption{Time of the Transition Between Phase 1 and Phase 2\tablenotemark{a}}
%\centering
%\begin{tabular}{l c}
%\hline
% Run  & Time (min)  \\
%\hline
%  $l1$  & 260   \\
%  $l2$  & 300   \\
%  $l3$  & 340   \\
%  $h1$  & 270   \\
%  $h2$  & 250   \\
%  $h3$  & 380   \\
%  $r1$  & 370   \\
%  $r2$  & 390   \\
%\hline
%\end{tabular}
%\tablenotetext{a}{Footnote text here.}
%\end{table}

% See below for how to make landscape/sideways figures or tables.

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

More Information and Advice:

%% ------------------------------------------------------------------------ %%
%
%  SECTION HEADS
%
%% ------------------------------------------------------------------------ %%

% Capitalize the first letter of each word (except for
% prepositions, conjunctions, and articles that are
% three or fewer letters).

% AGU follows standard outline style; therefore, there cannot be a section 1 without
% a section 2, or a section 2.3.1 without a section 2.3.2.
% Please make sure your section numbers are balanced.
% ---------------
% Level 1 head
%
% Use the \section{} command to identify level 1 heads;
% type the appropriate head wording between the curly
% brackets, as shown below.
%
%An example:
%\section{Level 1 Head: Introduction}
%
% ---------------
% Level 2 head
%
% Use the \subsection{} command to identify level 2 heads.
%An example:
%\subsection{Level 2 Head}
%
% ---------------
% Level 3 head
%
% Use the \subsubsection{} command to identify level 3 heads
%An example:
%\subsubsection{Level 3 Head}
%
%---------------
% Level 4 head
%
% Use the \subsubsubsection{} command to identify level 3 heads
% An example:
%\subsubsubsection{Level 4 Head} An example.
%
%% ------------------------------------------------------------------------ %%
%
%  IN-TEXT LISTS
%
%% ------------------------------------------------------------------------ %%
%
% Do not use bulleted lists; enumerated lists are okay.
% \begin{enumerate}
% \item
% \item
% \item
% \end{enumerate}
%
%% ------------------------------------------------------------------------ %%
%
%  EQUATIONS
%
%% ------------------------------------------------------------------------ %%

% Single-line equations are centered.
% Equation arrays will appear left-aligned.

Math coded inside display math mode \[ ...\]
 will not be numbered, e.g.,:
 \[ x^2=y^2 + z^2\]

 Math coded inside \begin{equation} and \end{equation} will
 be automatically numbered, e.g.,:
 \begin{equation}
 x^2=y^2 + z^2
 \end{equation}

% IF YOU HAVE MULTI-LINE EQUATIONS, PLEASE
% BREAK THE EQUATIONS INTO TWO OR MORE LINES
% OF SINGLE COLUMN WIDTH (20 pc, 8.3 cm)
% using double backslashes (\\).

% To create multiline equations, use the
% \begin{eqnarray} and \end{eqnarray} environment
% as demonstrated below.
\begin{eqnarray}
  x_{1} & = & (x - x_{0}) \cos \Theta \nonumber \\
        && + (y - y_{0}) \sin \Theta  \nonumber \\
  y_{1} & = & -(x - x_{0}) \sin \Theta \nonumber \\
        && + (y - y_{0}) \cos \Theta.
\end{eqnarray}

%If you don't want an equation number, use the star form:
%\begin{eqnarray*}...\end{eqnarray*}

% Break each line at a sign of operation
% (+, -, etc.) if possible, with the sign of operation
% on the new line.

% Indent second and subsequent lines to align with
% the first character following the equal sign on the
% first line.

% Use an \hspace{} command to insert horizontal space
% into your equation if necessary. Place an appropriate
% unit of measure between the curly braces, e.g.
% \hspace{1in}; you may have to experiment to achieve
% the correct amount of space.


%% ------------------------------------------------------------------------ %%
%
%  EQUATION NUMBERING: COUNTER
%
%% ------------------------------------------------------------------------ %%

% You may change equation numbering by resetting
% the equation counter or by explicitly numbering
% an equation.

% To explicitly number an equation, type \eqnum{}
% (with the desired number between the brackets)
% after the \begin{equation} or \begin{eqnarray}
% command.  The \eqnum{} command will affect only
% the equation it appears with; LaTeX will number
% any equations appearing later in the manuscript
% according to the equation counter.
%

% If you have a multiline equation that needs only
% one equation number, use a \nonumber command in
% front of the double backslashes (\\) as shown in
% the multiline equation above.

%% ------------------------------------------------------------------------ %%
%
%  LANDSCAPE/SIDEWAYS FIGURE AND TABLE EXAMPLES
%
%% ------------------------------------------------------------------------ %%
%
% For figures, add \usepackage{lscape} to the file and the landscape.sty style file
% to the paper folder.
%
% \begin{figure*}[p]
% \begin{landscapefigure*}
% Illustration here.
% \caption{caption here}
% \end{landscapefigure*}
% \end{figure*}
%
% For tables, add \usepackage{rotating} to the paper and add the rotating.sty file to the folder.
%
% AGU prefers the use of {sidewaystable} over {landscapetable} as it causes fewer problems.
%
% \begin{sidewaystable}
% \caption{}
% \begin{tabular}
% Table layout here.
% \end{tabular}
% \end{sidewaystable}
%
%